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ANALYSIS  OF  AN  ASSEMBLAGE  OF  DISCS 
EMPLOYING  INTERACTIVE  GRAPHICS 


PART  I:  INTRODUCTION 

Background 

1.  A  computer  modeling  -.echnique  termed  the  ’’distinct"  element 
method  was  introduced  by  Cundal 1  (1971).  Then,  several  years  later,  the 
technique  was  extended  to  analyze  blocky  rock  systems  (Cundall,  1974). 
This  type  of  analysis  permitted  large-scale  motions  of  the  individual 
rock  blocks.  This  work  was  further  extended  to  the  analysis  of  tunnel 
supports  (Voegele,  1979).  Thus,  the  distinct  element  method  has  been 
developed  to  accurately  model  many  features  of  the  behavior  of  jointed 
rock. 


Purpose 

2.  This  study  was  directed  at  developing  methods  for  analyzing 
the  behavior  of  systems  of  simple,  particulate  members.  The  basic  shape 
of  the  particulate  members  was  chosen  to  be  a  disc  (or  cylinder).  The 
distinct  element  method  is  oriented  more  toward  representing  the  kine¬ 
matics  of  a  system  of  particulate  members  instead  of  the  continuum 
aspects.  At  the  outset,  the  structure  is  viewed  as  a  collection  of 
individual  members  with  individual  properties,  rather  than  a  homogeneous 
system  for  which  the  material  property  effects  are  assumed  to  apply 
throughout . 

3.  In  addition  to  extending  the  distinct  element  method  to 
discs,  there  was  a  desire  to  use  computer  interactive  graphics  in  the 
analysis  so  that  the  user  could  be  relieved  of  laborious  data  preparation 
and  time-consuming  plotting  of  the  results.  This  report  describes  the 
development  of  a  computer  program  entitled  DISC  and  contains  discussions 
of  the  mathematical  formulations,  program  organization,  and  program 
operat ion. 
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PART  II:  MATHEMATICAL  CONCEPTS 

4.  The  computer  program  DISC  was  formulated  using  concepts  of 
the  "distinct  element”  method  proposed  by  Cundall  (1974).  Conceptually, 
the  distinct  element  method  computes  the  forces  and  displacements  at 
points  of  contact  between  mathematically  described  particles  (elements). 
The  contact  between  the  elements  are  mathematically  represented  by 
springs  and  dashpots.  Any  relative  displacement  between  two  elements  in 
contact  then  results  in  a  force  developing  within  the  springs  and  dash- 
pots.  The  vector  sum  of  all  such  forces  on  the  discs  in  turn  causes 
accelerations  of  the  discs.  If  these  accelerations  are  assumed  to  be 
constant  for  a  prescribed  interval  of  time,  the  associated  velocities 
and  displacements  may  be  computed  for  the  end  of  that  time  interval. 

These  new  displacements  may  then  be  used  to  determine  new  relative 
displacements  between  the  contacting  elements,  and  the  cycle  is  then 
repeated  over  and  over. 

5.  The  computer  program  DISC  is  capable  of  handling  element 
shapes  corresponding  to  discs  and  bars.  The  bar  elements  are  composed 
of  two  parallel  sides  of  arbitrary  length  terminated  on  both  ends  by 
semicircles.  Figure  1  shows  two  examples  of  two  elements  in  contact. 
Condition  a  presents  two  disc  elements  in  contact,  and  condition  b  a  bar 
and  a  disc  element  in  contact.  Mathematically,  the  equations  required 
to  describe  the  subsequent  motions  of  either  contact  pair  are  identical. 
However,  condition  b  is  somewhat  more  general  and  easier  to  visualize 
and  as  such  will  be  referred  to  in  the  following  discussion. 

6.  Figure  lb  shows  the  bar  element  i  ,  its  centroid  located  at 
x1,  y1  ,  and  the  disc  element  j  ,  its  centroid  located  at  x~* ,  y^  . 

The  coordinates  of  the  contact  point  are  at  x  ,  y  .  The  angle  to  the 
contact  plane  is  denoted  by  B  .  (For  the  conditions  shown  in  Figure  la, 
B  is  the  angle  of  the  plane  of  tangency.) 

7.  The  basic  key  equation  cf  the  distinct  element  method  relates 
the  incremental  displacement  of  the  contact  point  (into  one  of  the  two 
contacting  bodies)  to  the  relative  motion  of  the  elements  that  are  in 
contact.  The  following  equation  is  the  result  of  superposition,  i.e.. 
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Element  j 


Figure  1.  Basic  elements  and  definitions 


first  assume  element  i  is  fixed  and  determine  the  resulting  motion  of 
the  contact,  and  then  assume  element  j  is  fixed  and  examine  the  motion 
of  the  contact  point.  This  analysis  yields  (as  illustrated  in  Figure  2a) 


and 


where 


c 

x 

c 

y 

AO 


Au"  =  Au^  -  Au1  +  A(P(x°  -  x~* )  -  A0*(xC  -  x'S 

y  y  y 

(1) 

Au^  =  Au^_  -  Au^  -  AO*(yC  -  y^ )  +  A6^(yC  -  y*) 

the  vertical  component  of  the  incremental  contact 
d isplacement 

the  horizontal  component  of  the  incremental  contact 
d isplacement 

the  current  vertical  incremental  displacement  of  element  i  , 
etc.  (i.e.,  the  subscript  (x  or  y)  refers  to  the  component 
direction  and  the  superscript  (i  or  j)  refers  to  the  element) 

the  centroid  coordinates  of  element  i 

the  centroid  coordinates  of  element  j 

the  horizontal  coordinate  of  the  contact  point 

the  vertical  coordinate  of  the  contact  point 

the  incremental  rotation  of  the  element 


8.  The  global  angle  of  the  contact  plane  is  given  by  B  .  Thus, 
the  horizontal  and  vertical  relative  contact  displacements  may  be  re¬ 
solved  into  components  parallel  and  perpendicular  (shear  and  normal 
directions)  to  the  contact  plane.  Thus,  as  shown  in  Figure  2b 


and 


where 


c  c  c 

Au  =  Au  sin  B  +  Au  cos  B 
s  y  x 


c  c  c 

Au  =  Au  cos  B  -  Au  sin  B 
n  y  x 


(2) 


c 

Au  =  the  component  of  the  contact  point  displacement  parallel  to 
the  contact  plane 

Au^  =  the  normal  component 

9.  As  previously  mentioned,  the  contacts  are  represented  by 
springs  and  dashpots  (in  both  the  normal  and  shear  directions).  Thus, 
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(FC  +  nC)  cos  B  +  (f'L:  +  n1')  sin  B 
s  s  n  n 


y 

•j 

x 


Figure  2.  Contact  displacements  and  forces 


the  force  generated  within  the  springs  and  dashpots  beca'i1-".  ui  the 
incremental  normal  and  shear  displacements  is  given  by 


FC  =  l’C  -  :,uCV 
n  n  n  n 


l'C  =  i'C  +  AuCk 
s  s  s  s 


in*  K 
n  n 


I)1’  = 

s 


At/  K 
s  s 


(3) 


wiiere 


n 

F° 


1)C 

n 

nr 

s 

k 

n 

k 

s 

K 

n 
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t lie  normal  spring  force 

the  shear  spring  force 

the  normal  dashpot  force 

the  shear  dashpot  force 
the  normal  stiffness 
the  shear  stiffness 
the  normal  damping  constant 
the  shear  damping  constant 


The  symbol  =  means  "replaced  by";  that  is,  after  each  increment  of 
contact  point  displacement,  new  normal  and  shear  forces  are  recomputed 
as  functions  of  the  old  values.  The  dashpots  are  necessary  to  prevent 
the  two  contact  elements  from  vibrating  indefinitely.  To  be  strictly 
correct,  the  dashpot  forces  should  be  related  to  the  contact  point 
velocity;  however,  for  a  time  increment,  the  incremental  displacement  is 
proportional  to  the  velocity  (i.e.,  Au  =  vAt  ,  where  v  is  the  relative 
velocity  across  the  contact).  Figure  2c  shows  the  positive  directions 
of  the  normal  and  shear  forces.  Also,  note  that  a  positive  normal  force 
indicates  a  compressive  force. 

10.  Equation  '3  is  subject  to  modification  if  the  tensile 
"strength"  T  (a  negative  value)  is  exceeded.  Thus,  it  follows  that  if 
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element  i  ,  may  possess  a  number  of  contacts.  In  addition,  element  i 
may  be  subject  to  applied  forces  and  gravity  forces.  Thus,  the  total 
force  acting  on  element  i  is  given  by 


Fl  =  T.  F1  +  vertical  applied  load  -  gravity  load 
vsum  c  y 


F1  =  Y.  F1  +  horizontal  applied  load  (5) 

xsum  c  x 


=  >;  [fV  -  x^)  -  F^(yC  -  ySl  +  applied 

c  |_  y  x  J 


moment 


The  symbol  means  the  summation  for  all  contact  points  of  element  i  . 

The  equations  are  similar  for  element  j  .  The  term  M1  is  the  total 
n  J  sum 

moment  acting  on  element  i  .  Figure  3  illustrates  the  force  sums. 


ysum 


+  applied  moment 


Figure  3.  Force  sums  on  element  i 
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12.  After  computing  all  the  force  sums  acting  on  each  element, 
the  acceleration  of  each  element  may  be  computed  from  Newton's  law  of 
motion,  i.e.,  acceleration  =  unbalanced  force/mass,  and  angular  acceler¬ 
ation  =  unbalanced  moment/moment  of  inertia.  Thus,  the  component 
accelerations  a  ,  and  angular  acceleration  a  are  expressed  as 

i  _  ~i  ,  i 

a  =  F  /  m 

x  xsum 


a  =  F  /m 
y  ysum 


where 


i 

m 


i 

u 


sum 


the  mass  of  element  i 

the  mass  moment  of  inertia  of  element  i 


(6) 


The  velocities  (  to  is  the  angular  velocity)  are  then  computed  by 
multiplying  by  the  time  step  At  : 

i  .  i  ,  i, . 
v  =  v  +  a  At 

XXX 


v1  =  v1  +  a1 At 

y  y  y 


1  .  1  ,  1 A 

u>  =  ii)  +  i  At 


(7) 


A  second  numerical  integration  yields  the  incremental  displacement  for 
the  element  to  be 

Au1  =  v  At 
x  x 


Au 1  =  v  1 A  t 

y  y 


A  t 


.  l 


i  . 
c.L 


(«) 


The  quantities  above  are  sought  in  order  to  recycle  the  calculations. 
These  quantities  are  introduced  back  into  Equal  ion  1,  and  the  whole 
process  is  repeated. 
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13.  Finally,  the  total  displacement  of  the  element  centroids  are 
found  from 

i  i  i  ,  ,  i 
u  =  u  +  Au 

XX  X 

u1  =  u1  +  Au1  (9) 

y  y  y 

o1  =  o1  +  A<’1 

14.  At  the  outset  of  each  problem,  the  incremental  displacements 
are,  of  course,  zero.  Thus,  Equations  1  through  4  yield  nothing,  i.e., 
there  are,  up  to  now,  no  contact  forces.  In  Equation  5,  a  force  sum 
will  result  due  to  applied  and/or  gravity  forces.  Thus,  an  acceleration 
of  an  element  will  occur  that  in  turn  can  yield  incremental  displacements. 
Equations  1  through  4  are  then  used  to  compute  the  new  contact  forces. 

15.  As  noted  from  the  discussion  above,  the  mathematics  involved 
in  program  DISC  are  quite  simple,  since  only  the  concepts  of  a  damped 
oscillator  and  Newton's  law  of  motion  are  needed.  The  real  difficulty 
in  performing  this  type  of  analysis  is  the  development  of  the  logic  to 
detect  and  to  provide  computer  memory  for  the  many  contacts  that  may 
result.  Thus,  the  problem  becomes  one  of  geometry  and  bookkeeping.  The 
details  of  the  program  are  presented  in  Part  Ill. 
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PART  111:  PROGRAM  ORGANIZATION 


General  Overview 


16.  rin-  computer  program  DISC  is  a  time-sharing  FORTRAN  program 
lliat  consists  of  a  main  program  and  sixteen  subroutines.  Physically, 

I)  i  SO.  is  composed  ol  approximately  1200  FORTRAN  statements.  The  program 
is  written  lo  use  interactive  graphics  and  is  operated  by  a  Tektronix 
4014  (or  4010)  termin.il.  With  the  except  ion  of  occasionally  tvping  in  a 
few  numbers  to  rcdei  ine  certain  parameters,  the  commands  primarily 
consist  ot  single  keystrokes.  flic  output  is  primarilv  graphical,  i.e., 
pictures  ol  the  present  local  ions  ol  the  elements  are  drawn  on  the 
screen.  While  it  is  possible  to  obtain  printed  information,  the  normal 
mode  of  operation  is  visual. 

17.  As  discussed  in  Part  II,  the  actual  calculations  lor  deter¬ 
mining  the  mot  ion  ol  th«-  disc  and  liar  elements  are  quite  s  t  ra  ight  f  "rwa  rd . 
Knowing  the  location  ol  a  contact  poinL  and  the  previous  relative  d i s- 
placenient  s  of  those  elements  in  contact,  I  lie  lorces  between  the  elements 
can  be  computed  and  then  integrated  over  time  to  obtain  new  displacements 
neeessarc  tor  the  next  evcle.  This  process  is  carried  out  lot  til’, 
contacts  at  each  time  ivcle. 

IS.  Suppose  though,  that  a  system  of  100  disc  elements,  packed 
closely  together,  is  to  he  onsidered.  II  the  discs  arc  of  the  same 
size,  then  there  could,  at  most,  be  six  contacts  per  disc;  and  since 
some  must  he  on  or  near  the  edge,  I  lie  total  number  of  contact  points 
must  be  less  than  POO.  However,  in  older  to  chick  for  coiilui  ts  between 
the  two  arbitrary  discs,  it  would  (bv  brute  force)  he  ntu  essarv  to  check 
for  the  possibility  of  a  contact  lie  tween  each  and  every  disc.  I'utllur- 
more,  an  inordinate  amount  ol  searching  would,  be  required.  !'  n  100 
discs,  the  number  of  litlerent  possible  combinations  i  is  civcn  bv 


( 


n! 

2  (  n  -  2)  ! 


100! 
2(()8! ) 


49  60 
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For  300  discs,  the  number  of  searches  would  increase  to  124,750.  In 
addition,  computer  memory  space  would  be  required  at  double  the  amounts 
above  to  store  the  normal  and  shear  contact  forces. 

19.  Since  most  ot  this  memory  would  be  blank  and  most  of  the 
searches  fruitless  (where  contacts  were  not  found),  a  great  deal  of 
effort  was  put  into  formulating  schemes  and  methods  to  facilitate  the 
program's  efficiency.  Indeed,  the  usability  of  the  distinct  element 
method  is  predicated  on  elficient  programming  techniques.  Without 
incorporating  these  schemes  it  would  not  be  possible  to  solve  even  small 
problems  on  large  computers. 

20.  The  various  subroutines  that  make  up  the  program  DISC  are 
related  to  each  other  primarily  through  named  common  storage.  The  items 
in  common  storage  are  as  follows: 

N  '■  Number  of  disc  and  bar  elements 
SKN  =  Spring  stiffness  at  contact 
DKN  =  Damping  constant  of  contacts 
TIMK  *  lime  since  beginning  of  run 
DT  =  Tiini’  increment  between  steps 

DF  =  Damping  factor  (a  multiplier  for  the  damping  constant) 

FXSliM(M)  =  Sum  of  forces  iti  x-direction  on  element  M 

FYSIIM(M)  =  Sum  of  forces  in  y-diroct.  ion  on  eleme,.'  M 

FMSUM(M)  =  Sum  of  moments  acting  on  element  M 

FX(.M)  Applied  x- force  on  element  M 

FY(M)  -  Applied  y-force  on  element  M 

F>1(M)  =  Applied  moment  on  element  M 

X(M)  =  Current  x-coard inate  of  element  M 

Y  (M)  =  Current  y-eoord inate  of  element  M 

T ( M )  =  Current  angular  rotation  of  element  M 

K(M)  =  Radius  (or  thickness,  if  bar  element)  of  element  M 

S(M)  =  One  half  of  the  separation  between  the  circular  ends 
of  bar  element  M 

W(M)  -  Weight  of  element  M 

XM(M)  =  Mass  of  element  M 

UX(M)  -  Current  x-d i sp 1 acement  of  element  M 


LIY(M)  = 
UT(M)  = 
VX (M)  = 
VY  (M)  = 
1FIX(M)  = 


wc;  ri'(M)  = 

FRC(M)  = 

TKN(M)  - 
AA(M  x  24)  = 


NA(M)  = 


1,1.  (M  x  ft)  - 
The  common  S(.orngo 


Current  y-d isp 1 acemenL  of  element  M 

Current  angular  d isp 1 ae ement  of  element  M 

Current  x-veloeity  of  element  M 

Current  y-veloeitv  of  element  M 

Fixity  code  for  element  M  .  The  values  are: 

1  =  Fixed  in  x-directiou 

2  =  Fixed  in  y-direetion 

1  =  Fixed  in  x-  and  y-direction 

4  -r  Fixed  in  angular  direction 

5  -  Fixed  in  x  and  angular  directions 

b  =  Fixed  in  v  and  angular  directions 

/  -  Fixed  in  all  directions 

8  =  A  constant  preset  ihed  velocity  applies 

Weight  (actor  (a  multiplier  for  the  initially 
assigned  weights) 

Angle  corresponding  to  coefficient  of  friction, 
i . e.  tan( FKC)  =  n 

TonsiLc  strengln  at  contact  for  element  M 

AA  is  tiie  contact  list.  See  discussion  of 
subroutine  MOTION 

Index  used  to  find  sought  quantities  in  the 
eontaet  list 

1’ossible  contact  search  list.  See  discussion  of 
subrout  ine  FI.  i  S  !' 

requites  approxim.it  el  v  S3  N  storage  locations. 


Main  1’rogram 

21.  The  main  program  initiates  the  graphics  software  and  calls 
subroutine  C1DYUP.  Control  is  never  returned  to  the  main  program.  The 
graphics  subroutines  (West  Point  Military  Academy,  1975)  are  peculiar  to 
the  WF.S  computing  system  and  also  the  Office  of  Personnel  Management's 
Honeywell  system  at  Macon,  Ocorgia.  The  graphics  subroutines  employed 
in  DISC  include: 
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U START 

UHKLL 

UCRCLE 

USET 

U ERASE 

Ul’EN 

UWA  1 T 

UPRINT 

UMOVE 

l.'K  FAD 

UPRNT1 

UPSET 

UCKIN 

UARC 

UDAREA 

l!W  1 NDO 

U ROTATE 

UOUTLN 

Subroutine  OIUYUP 


22.  Subroutine  OIDYUP  is  the  driving  subroutine.  Jts  primary 
function  is  to  call  subroutine  INPUT  and  to  cause  cycling  through  sub¬ 
routine  OUTPUT.  This  subroutine  is  the  only  one  that  calls  the  main 
calculation  subroutine  MOTION.  The  normal  sequence  of  operations 
involves: 

a.  OIDYUP  calls  INPUT.  The  user  interactively  describes 
the  problem.  Alternatively,  the  user  may  select  a 
problem  previously  stored  by  subroutine  SAVE. 

b.  A  call  is  made  to  PI, 1ST,  which  prepares  a  list  of  those 
elements  tor  which  a  contact  is  possible  within  the  next 
50  cycles. 

c.  A  call  is  made  to  MOTION,  which  calculates  the  movements 
of  and  forces  on  the  elements.  Fifty  cycles  are  per¬ 
formed  in  MOTION.  At  the  end  of  each  50  cycles,  the  new 
positions  of  the  elements  are  graphically  plotted  by 
OUTPUT.  Step  b  is  again  performed. 

d.  Every  so  often  (user  selected,  default  value  =  250  cycles) 
a  call  is  made  to  OUTPUT.  The  program  operation  is 
halted,  and  the  user  must  inform  the  system  what  action 

is  to  be  taken  next  (i.e.,  continue  running,  draw  all 
elements,  change  parameters,  etc.). 


Subroutine  INPUT 


23.  Subroutine  INPU'I  is  used  to  describe  the  problem  to  be 
analyzed  (see  user's  guides  (Figures  9  and  10))  and  also  to  initialize 
common  data  storage  dependent  upon  its  calling  subroutine.  It  may  be 
called  by  OIDYUP  (to  start  a  new  problem)  or  by  OUTPUT  (to  modify  an 
existing  problem).  The  following  subroutines  may  be  called  by  INPIT: 
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CIRCLE 
I  OOP 
CEN 


SCREEN 
C1DYUP 
OR  ID 
SAVE 


Subroutine  OUTPUT 

2-4.  A  variety  of  functions  can  be  performed  by  subroutine  OUTPUT 
(see  user's  guides  (Figures  9  and  10)).  After  the  problem  geometry  is 
described  in  INPUT,  subroutine  OUTPUT  takes  over  the  system  operation. 

Its  function  is  to  set  boundary  conditions  (fix  certain  elements),  to 
allow  for  inputting  parameters  (weights,  damping  factors,  frictional 
properties,  etc.),  and  to  handle  both  graphical  and  formatted  output. 
Between  the  phases  when  calculations  of  motion  are  being  performed 
between  CIDYUP  and  MOTION,  the  program  halts  occasionally  (user  defined, 
default  value  =  250  cycles)  in  OUTPUT.  During  this  halt,  the  user  may 
redefine  boundary  conditions,  parameters,  etc.,  plot  the  data,  or  request 
another  calculation  phase. 

25.  During  any  program  halt,  the  user  may  also  instruct  the 
program  to  store  a  copy  of  the  existing  data  geometry  (performed  by 
subroutine  SAVE)-  Ibis  stored  data  may  then  be  recalled  by  a  direct 
call  to  INPUT.  Subroutine  OUTPUT  is  also  used  to  delete  (erase)  certain 
elements;  however,  elements  may  not  be:  added  in  OUTPUT.  At  the  outset 
of  a  problem,  a  great  deal  of  geometrical  editing  may  be  accomplished  by 
juggling  the  program  back  and  lortb  between  OUTPUT  and  INPUT.  Conse¬ 
quently,  the  user  is  able  to  reposition  elements,  add  elements,  delete 
elements,  change  boundary  conditions,  etc. 

2b.  The  data  storage  (SAVE)  teature  is  very  handy  for  situations 
in  which  it  is  desired  to  solve  a  suite  of  problems  with  the  same  geometry 
but  difterent  parameter' .  Alter  one  problem  is  solved,  the  stored  data 
may  he  recalled,  t  lie  parameters  changed,  ami  the  next  problem  solved. 

27.  Subroutine  OUTPUT'  mav  be  called  only  by  subroutine  OIDYUP. 

The  following  subroutines  may  be  called  by  OUTPUT. 

CIDYI’P  (  I  RULE  DAMP 

INPUi  VECTOR  I’Ll  TOR 

LOOP  WE  I  CUT  OK  1 D 

SCREEN  INTERVAL  SAVE 


28.  Subroutine  OUTPUT  is  also  used  to  compute  the  contact  stiff¬ 
ness  k  and  time-step  increment  At  .  (Subroutine  INPUT  tentatively 
sets  k  and  At  in  the  same  fashion  as  OUTPUT.  However,  OUTPUT 
redefines  a  stable  time  step  dependent  upon  the  boundary  conditions.) 

29.  The  equation  of  motion  for  an  undamped  linear  oscillator  of 
stiffness  k  and  mass  m  is  given  by 

.2 

m  — —  +  ku  =  0 
dt 

where  u  is  the  displacement.  In  a  finice  difference  notation,  the 
2  2 

derivative  d  u/dt  may  be  expressed  as 


dju 

dt 


u  ,  ,  -  2u  +  u 
t+1  t  t-1 


At 


where  ut+1  >  ut  >  and  ut_|  represent  the  displacements  at  times 
t+1  ,  t  ,  and  t-1  ,  respectively.  Substitution  into  the  equation  of 
motion  results  in 

ut+l  +  (m  At2  "  2)  Ut  +  Vl  =  0 

The  equality  of  the  finite  difference  equation  above  can  hold  only  if 
the  quantity  |  At2  -  2j  is  negative.  Thus,  the  stable  time  step  At 
is  given  by 


At  < 


(11) 


Greater  values  for  At  will  lead  to  exponentially  increasing  u  . 

30.  It  can  be  shown  that  the  formulation  for  DISC  (when  consid¬ 
ering  a  single  normal  contact)  yields  the  same  finite  difference  equa¬ 
tion  given  above.  That  is,  given  (for  some  time  step  t)  the  values 
of  ut  ,  Au^  ,  vf  ,  and  F  in  the  normal  directions.  Equation  3 
becomes 


F 


t+1 


t 

-  kAu  =  ->;  kAu  =  -ku 
C  t=0  C 
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llquat  ions  b  and  7  may  bo  expre  ssed  as 


!•  ,  ,  At  ,  , 

t+1  kuAt 

v  ,  =  v  +  -  -  v  -  —  — 

t+1  t  in  t  m 


From  liquation  8  and  by  substitution  nt  tlio  preceding  equation 

2 

ku.it  ~ 

A u  ,  =  v  , .At  =  v  At  -  - 

t  +  1  t+1  L  ill 


Then,  making  the  substitution 


v  . .  t  -  ai  =  u  -  u 
t  t  r  i  -i 


resu Its  in 


ku  .  t 
L 


'l't+l  “t  ''(-I 


Finally,  Equation  9  may  be  written  as 


U  ,  ,  =11  i 

t  +  i  t 


a  + 
t 


(  ,!l  +  I  ) 


-  u  , 
t  t-1 


ku  \L 

t 


(- 


Vt 


u  -  u 
t  t-1 


Or  by  rearranging  terms 


u  +  (-•  At  —  \  u  -  u , 

t+j.  \tn  !  t  t-1 


The  displacement  at  time  t+i  can  be  computed  i!  the  two  previous 
displacements  at  times  ;  and  t-1  are  known. 

if.  The  procedure  outlined  for  the  computat  ions  of  a  stable  time 
step  At  considered  on  I v  one  contact.  In  actuality,  many  contacts  can 
be  made  for  a  single  element ;  thus,  the  useful  time  stop  must  be  chosen 
smaller  than  the  single  contact  time  step.  That  is,  most  problems 
consist  of  many  degrees  of  freedom. 


I  9 


32.  Thus,  a  stable  At  will  increase  as  the  mass  increases  and 
will  decrease  as  the  stiffness  increases.  The  program  will,  of  course, 
calculate  faster  (i.e.,  the  cost  will  go  down)  as  the  time  step  increases. 
Alternatively,  if  a  At  is  chosen,  then  a  stable  k.  could  be  chosen 

for  a  given  m  .  In  general,  any  two  of  At  ,  k  ,  and  m  can  be 
chosen,  and  the  other  computed  to  ensure  stability  through  the  preceding 
equation.  The  scheme  incorporated  into  program  DISC  for  computing  a 
stable  At  is  described  in  the  following  paragraphs.  Any  subsequent 
user  of  DISC  should  not  necessarily  feel  bound  to  this  scheme. 

33.  At  the  outset  of  writing  DISC,  it  was  decided  that  a  disc 
element  with  a  radius  of  20  screen  units  would  be  the  "standard"  disc 
element.  The  Tektronix  screen  is  1023  screen  units  wide  and  780  screen 
units  high.  For  this  standard  element,  a  weight  of  100.0  units  is 
assigned.  To  determine  the  weight  of  other  disc  and  bar  elements,  the 
area  of  the  element  is  computed,  divided  by  the  standard  area,  and 
multiplied  by  100.0.  Thus,  all  forces  output  by  the  program  are  in 
terms  of  the  100.0  units  assigned  to  the  standard  disc,  i.e.,  the  forces 
are  normalized  and,  as  such,  are  nond imensional .  The  element  mass  m 

e 

is  determined  from 


ne  =  ue^  U-2) 

where 

=  the  normalized  weight  of  the  element 

g  =  the  acceleration  due  to  gravity 

The  value  g  is  an  input  parameter  (default  value  =  32.2).  Thus,  if 

the  default  value  of  32.2  is  used,  it  implies  an  acceleration  due  to 

gravity  of  32.2  screen  units/unit  of  time/unit  of  time.  If  time  is 

reckoned  in  seconds,  this  would  normally  suggest  that  one  screen  unit 

2 

equals  one  foot  (i.e.,  g  =  32.2  ft/sec  ).  Inputting  g  as  9.8  (with 

time  reckoned  in  seconds)  suggests  that  one  screen  unit  equals  one  metre 

2 

(i.e.,  g  -  9.8  m/sec  ).  It  is  through  this  input  parameter  g  ,  that  a 
scale  can  be  assigned  to  the  problem. 
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34.  I'lie  i'  out  net  spring  si  i !  I  ness  k  (cu  r  n  ■  ill  1  y  set  equal  in 
normal  and  shear  d i ree t ions)  can  be  assigned  in  a  variety  of  ways.  If 
the  stillness  were  physically  known,  i!  <  ould  he  din-el  ly  assigned.  In 
most  quasi— sta Lie  problems,  however,  the  purpose  ot  the  spring  stiffness 
is  to  prevent  the  elements  from  peneLrating  Loo  far  into  each  other. 

This  requirement  implies  choosing  a  very  large  k  ;  however,  choosing  a 
large  k  results  in  a  small  t  .  For  program  DISC,  the  following 
equation  is  used  for  selecting  k  : 


k  -  v  yJx  +  i  (ii) 

max 

where 

N  =  the  number  of  element  a  i a  the  system 

w  =  the  weight  ol  the  largest  clement  (actually  the  weight 
ot  the  largest  element  not  restrained  from  all  movement 
or  "fixed") 

31.  file  time  step  t  is  then  computed  as 


M  =  0.05\/m  .  /k  (14) 

»  m  ui 

where  r.i  .  is  the  mass  of  the  sma  i  lest  element  not  restrained  from  all 
min 

movement.  The  fac  tor  0.03  is  used  to  ensure  stability  since  more  than 
one  contact  per  element  is  possible.  In  addition,  subroutine  OUTPUT 
computes  element  areas  and  moments  of  inertia  and  initializes  many  data 
1  ists . 


Subroutine  1'hlST 


36.  Subroutine  PI. 1ST  is  used  to  prepare  a  data  list  of  possible 
contacts  (between  elements).  This  subroutine  is  ui cessed  bet  ore  each 
call  to  MOTION  (i.c.,  every  30  cyci.es).  Suppose  tin-  >•  i  tuaf  ion  as  shown 
in  Figure  4  exists.  To  del  ermine  which  elements  would  go  into  the 
possible  contact  list,  a  check  is  t  it's l  made  on  Liu-  distance  tL  . 
between  each  element,  i.e.. 
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V  =  SOAtfv1  -  ) 

X  XX 


V.  =  50At  (v 1  -  v J  ) 

j  y  y 


Now,  if  d..  '  de.  (i.e.,  che  original  distance  apart  is  less  than  the 

extrapolated  distance  apart),  then  the  element  pair  is  moving  apart  and 

G 

the  contact  is  not  possible.  However,  if  d..  >  d..  ,  the  element  pair 

ij  iJ 

is  closing  and  this  pair  is  entered  into  the  possible  contact  list 

provided  that  d? .  <  RNG..  . 

ij  1J 

38.  At  first  sight,  it  would  appear  that  the  last  test  is  all 

that  is  necessary  to  accept  or  reject  the  contact  candidate.  In  fact, 
it  is;  however,  to  speed  up  computat ions ,  it  is  more  efficient  to  check 
for  possible  contacts  in  the  order  described.  As  soon  as  a  contact 
candidate  is  accepted  or  rejected,  additional  computations  (involving 
the  velocities)  are  no  longer  necessary.  The  idea  is  to  accept  or 
reject  with  as  few  computations  as  possible.  It  must  be  mentioned  that 
PLIST  forms  only  a  possible  contact  list.  This  list  is  used  in  subrou¬ 
tine  MOTION,  where  a  closer  examination  is  made  to  determine  if  contact 

is  actually  made.  For  example,  as  shown  in  Figure  4a,  the  bar  element  2 
obviously  does  not  make  contact  with  disc  element  1.  However,  it  will 
be  included  in  the  possible  contact  list.  Also,  elements  1  and  3  may, 

if  they  are  moving  together  rapidly  enough,  be  put  in  the  possible 

contact  list,  but  they  may  or  may  not  physically  come  into  contact 
during  the  next  50  cycles. 

39.  Searching  for  possible  contacts  is  done  in  the  following 
manner.  Starting  with  element  i  =  1  and  continuing  incrementally  to 

i  =  n  (the  total  number  of  elements),  distance  checks  are  made  with  all 
elements  j  with  j  ranging  from  i  +  1  to  n  .  Whenever  a  possible 
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contact  with  element  i  is  found,  the  value  /  (the  number  oi  the 
element  contacting  i)  is  imbedded  in  the  contact  list.  Alter  all 
possible  contacts  with  element.  i  have  been  found,  the  value  i  is 
imbedded  in  the  list  to  signal  the  end  of  all  contacts  with  element  i  . 
For  example,  if  it  is  assumed  that  elements  1  and  i,  shown  in  Figure  4a, 
cannot  make  contact  during  the  next  it)  cedes,  tin  possible  contact  list 
will  be 

2  —  4-1-2-'*-  i  -  '* 

Thus,  element  I  makes  possible  contacts  with  2  and  ,  element  2  makes  no 
new  contacts  (the  1-2  contact  has  already  been  found),  dement  3  makes 
contact  with  4,  and  of  course,  element  4  could  not  possibly  make  any  as 
yet  undetected  contacts.  The  underlined  numbers,  1,  2,  3,  and  4,  are 
the  ascending  values  of  i  and  signify  the  end  ol  each  segment  of  the 
list.  Now,  if  it  is  assumed  that  elements  1  and  3  could  make  contact, 
the  list  would  be 

2  -  3-  4-1-  2  -  4-  3-  4 

For  illustrative  purposes,  consider  the  arrangement  of  elements  shown  in 
Figure  4b.  The  possible  contact  list  for  the  system  would  be 

2  —  3—  4  —  5  —  f>  —  I  —  2  —  4—  7  —  _J—  5—  7  —  45  —  4  — 

b  -  8  -  9  -  5  -  9  -  h  -  8  -  7  -  9  -  8  -  9  -  11  -  12  -  10  - 

12  -  11  -  12 

or  a  total  of  19  contacts.  Notice  that  the  numbering  is  always  forward. 
There  is  no  reason  to  search  backwards  as  those  contacts  will  already 
have  been  found.  The  list  segments  are,  in  general,  longer  for  tire 
lower  numbered  elements.  Within  subroutine  PLIST,  the  possible  contact 
list  is  entitled  Id.  (Mi  1ST).  1  he  reserved  length  of  LI.  Is  six  times  the 

number  of  elements,  which  should  be  adequate  tor  most  problems  since 
this  length  will  handle  an  average  of  six  contacts  per  element.  There 
is  the  possibility  that  the  mentor v  reserved  could  be  exceeded  if  a  large 
number  of  bar  elements  are  present  in  the  problem  to  lie  analyzed. 
Seb'Mutine  FUST  is  called  only  bv  c.'lDVhP.  No  other  subroutines  are 
called  by  PI,  1ST. 


25 


Subroutine  MOTION 


40.  Subroutine  MOTTON  performs  three  basic  functions:  (a)  it 
detects  element  contacts,  (b)  it  computes  the  forces  acting  on  each 
element,  and  (c)  it  computes  the  subsequent  displacement  of  the  element 
due  to  those  forces.  Each  entry  to  MOTION  results  in  50  time  step 
iterations.  Almost  all  of  the  processor  time  used  during  the  operation 
of  DISC  is  used  in  this  subroutine. 

41.  Figure  5  is  the  flowchart  for  subroutine  MOTION.  As  pre¬ 
viously  mentioned  in  the  discussion  of  subroutine  PLIST,  a  data  list  of 
possible  contacts  is  created  by  PLIST  before  each  call  to  MOTION.  This 
list  contains  the  elements  that  are  presently  touching  or  could  be  close 
enough  to  touch  during  the  next  50  time  cycles  (i.e.,  the  duration  for 
each  call  to  MOTION).  While  within  MOTION,  only  those  contact  pairs 
contained  in  the  possible  contact  list  are  considered  as  potential  real 
contacts.  It  is  not  difficult  to  visualize  that  this  list  will  be  quite 
small  compared  to  a  list  unrestricted  by  distance  considerations. 

Indeed  it  is  likely  (at  least  for  disc  elements)  that  the  possible 
contact  list  is  a  close  approximation  to  the  actual  contacts. 

Conditions  for  contacts 

42.  Figure  6  shows  the  various  conditions  for  contacts  between 
disc  and  bar  elements.  When  describing  element  contacts,  the  convention 
is  to  state  that  element  i  contacts  element  j  ,  where  j  is  always 
greater  than  i  .  Searches  for  contacts  are  always  forward,  i.e.,  the 
lower  numbered  element  is  said  to  contact  the  higher  numbered  element, 
and  never  vice  versa. 

43.  In  general,  four  situations  must  be  considered  to  detect 
jreal  element  contacts  from  the  possible  contact  candidates  i  and  j 
provided  by  PLIST. 

a.  Both  elements  i  and  j  are  disc  elements  (Figure  6a); 
only  one  contact  can  result. 

b.  Element  i  is  a  bar  element  and  j  is  a  disc  element 

(Figure  6b);  only  one  contact  can  result. 

c.  Element  j  is  a  bar  element  and  i  is  a  disc  element 

(Figure  6c);  only  one  contact  can  result. 
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Figure  6.  Detection  of  contacts 


d.  Both  elements  1  and  j  are  bar  elements  (Figures  6d  , 
e,  f,  g,  and  h) ;  one,  two,  or  four  contacts  can  result. 

44.  For  situation  a,  it  is  simply  necessary  to  determine  if  the 
disc  elements  are  within  the  proximity  criteria  (to  be  discussed  subse¬ 
quently)  in  order  to  be  real  contacts. 

45.  For  situation  b,  determine  if  disc  element  j  is  within  the 
proximity  criteria  of  the  periphery  of  bar  element  i  .  If  so,  a  real 
contact  is  detected. 

46.  For  situation  c_,  determine  if  disc  element  i  is  within  the 
proximity  criteria  of  the  periphery  of  bar  element  j  .  If  so,  a  real 
contact  is  detected. 

47.  For  situation  d_,  it  is  necessary  to  restrict  the  location  of 
the  points  of  contacts  whenever  two  bar  elements  are  parallel.  For  two 
parallel  elements,  the  contact  points  are  defined  to  be  at  the  locations 
depicted  in  Figures  6d ,  e,  f,  g,  and  h,  i.e.,  where  circular  ends  of  the 
elements  become  tangent  to  the  central,  straight  portions.  All  of  the 
contacts  may  be  detected  by  considering  two  conditions:  (1)  element  i 
to  be  a  bar  element  and  element  j  to  be  composed  of  two  disc  elements 
(the  discs  being  located  at  either  end  of  bar  j  ),  and  (2)  element  j 

to  be  the  bar  element  and  element  i  to  be  composed  of  two  disc  elements. 

48.  Thus,  as  noted  in  Figure  6d ,  the  necessary  two  contacts  will 
be  detected  if  i  is  considered  to  be  the  bar  (condition  (1)  above)  and 
j  to  be  the  two  discs  since  both  circular  ends  of  bar  j  touch  bar  i  . 
However,  no  contacts  would  be  detected  if  j  is  the  bar  (Figure  6e, 
condition  (2)  above)  since  neither  end  of  bar  i  makes  contact  with  bar 
j  .  Figures  6f  and  6g  indicate  that  one  contact  is  detected  by  choosing 
i  as  the  bar  and  the  other  found  by  choosing  j  as  the  bar.  Figure  6h 
represents  the  special  case  of  two  bars  of  equal  length  parallel  to  one 
another.  In  this  case,  four  contacts  (two  coincident  pairs)  are  detected; 
two  from  choosing  i  as  the  bar  and  two  from  ciioosing  j  as  the  bar. 

49.  For  the  situation  in  which  two  bar  elements  are  not  parallel, 
then  only  one  contact  between  them  is  possible.  However,  the  scheme 
just  outlined  for  detecting  contacts  between  bar  elements  is  also  used 
when  examining  skewed  bar  elements.  The  proximity  criteria  will  reject 
any  potential  contacts  not  sufficiently  close. 


50.  It  must  be  mentioned  that  no  matter  whether  element  i  or 

j  or  both  are  considered  as  the  bar  element  (and  the  other  as  composed 
of  two  discs)  when  seeking  contacts,  the  particulars  about  that  contact 
are  always  associated  with  the  lower  numbered  element,  i.e.,  element  i 
This  procedure  is  in  keeping  with  the  "going  forward"  rule. 

Proximity  crite r  i  a 

51.  The  criteria  for  accepting  or  rejecting  a  possible  contact 
are  as  follows: 


a.  If  the  peripheries  of  the  two  elements  are  within  two 
screen  units  of  each  other  at  the  subject  contact 
point  (recall  that  the  "standard"  disc  element  radius 
is  20  screen  units),  the  contact  is  accepted  as  a 
cand idate  contact.  Otherwise,  it  is  rejected  for  this 
iteration  cycle. 

b.  After  a  possible  contact  is  elevated  to  cand Idate 
contact  status,  two  more  tests  are  necessary  to  finally 
accept  the  contact  as  a  real  contact.  If  this 
candidate  contact  is  present  in  the  contact  list 
generated  at  the  last  time  step  cycle,  it  is  accepted 
as  a  real  contact  and  put  in  the  new  real  contact 
list.  If  the  cand idate  contact  is  not  present  in 

the  previous  contact  list,  then  the  two  elements 
must  be  within  0.1  screen  units  of  each  other  in 
order  to  be  placed  in  the  new  real  contact  list. 
Otherwise,  it  is  rejected  for  this  iteration  cycle. 

Contact  list 

52.  The  contact  list  is  a  vector  array  in  COMMON  computer  stor¬ 
age.  This  list  contains  the  information  about  each  contact  that  is 
necessary  for  further  processing  through  subroutine  MOTION.  The  contact 
array  is  defined  by  the  named  COMMON  storage  statement: 

C0MM0N/ST0RE/AA(M  x  24),  NA(M) 

Since  M  is  the  maximum  number  of  elements  (N  is  the  actual  number  of 
elements)  that  may  be  considered  for  a  problem,  the  size  of  AA  is  re¬ 
served  to  be  at  24  times  M  .  This  size  for  AA  should  be  adequate  for 
most  problems  since  only  actual  contacts  are  stored  at  any  given  step  in 
the  calculations.  The  contents  of  AA  are  constantly  changing  as  the 
calculations  proceed.  The  array  AA  should  be  thought  of  as  divided  into 
two  segments:  the  first  segment  containing  the  contact  information  for 
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t  lu'  just  completed  iteration  (tin-  previous  contai  L  list),  and  the  second 
.segment  containing  the  intoimation  for  the  current  contact  list  (or  the 
now  contact  list).  As  contacts  are  detei  ted,  the  contact  information  is 
embedded  consecutively  within  the  contact  list.  Since  a  given  element 
can  contact  numerous  other  elements  at  any  given  time,  the*  consecutive 
nature  ot  the*  list  is  important  to  conserve*  memory.  flu*  array  NA(i) 
given  above  is  used  to  point  to  the  location  within  AA  where  the  con¬ 
tacts  for  a  given  element  i  begin. 

5  1.  Before  describing  the  actual  structure  of  the  contact  list, 
a  discussion  of  the  particular  quantities  going  into  the  contact  list  is 
beneficial.  For  ex. unpie,  suppose  a  contact  list  is  being  built  and  the 
next  available  location  within  the  array  AA  is  AA(K) ,  where  K  is  Lhe 
index  within  common  storage.  Further,  suppose  that  a  contact  has  just 
been  detected  between  elements  i  and  j  and  accepted  as  a  real  con¬ 
tact.  Also,  suppose  that  this  same  contact  was  present  in  the  previous 
contact  list  starting  at  some  locat  ion  AA(L).  Thus,  the  following 
information  is  put  into  the  new  contact  list  starting  at  AA(K): 

AA ( K )  =  j  +  v  /J 0000 

c 

AA  (K  +  1)  =  II)  +  v  / 10000 
c 

AA(K  +  2)  =  AA(h  +  2),  the  previous  normal  contact  force 

AA(K  +  J)  =  AA(I.  +  1),  the  previous  shear  contact  force 

AA(K  +  4)  =  B,  the  angle  from  the  horizontal  to  the  contact 
p  1  ane 

The  symbols  x  and  y  are  the  coordinates  of  the  point  of  contact, 
and  ID  is  an  identifier  indicating  the  fashion  in  which  the  contact 
was  detected  (i.e.,  for  all  cases  except  two  bar  elements  in  contact, 

ID  =  1  ;  for  two  bar  elements,  II)  takes  on  values  of  1  to  4  depending 
in  which  sequence  the  contacts  were  found  during  the  search).  In  effect, 
locations  AA(K)  and  AA(K  +  1)  each  contain  two  pieces  of  information. 

That  is,  AA(K)  contains  as  its  whole  part  the  element  number  j  (con¬ 
tacting  element  i)  and  as  its  fractional  part  x  .  A  similar  situation 

holds  for  AA(K  +  1). 

54.  Strictly  speaking,  only  the  local  ions  AA(K  +  2)  and  AA(K  +  1) 
contain  information  (the  normal  and  shear  forces  on  the  contact)  that 
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must,  be  retained  between  iterations  in  order  to  accomplish  the  calcula¬ 
tions  lor  displacement.  The  need  ior  storing  the  variables  j  and  ID 
is  to  a  itl  in  determining  it  the  contact  now  being  considered  is  present 
in  the  previous  contact  list.  The  variables  x  ,  v  ,  and  B  are 
retained  in  the  contact  list  as  an  aid  lor  later  outputting  of  the 
resul t  s . 

55.  It  it  is  now  supposed  that  the  contact  being  processed  was 
not  found  present  in  the  previous  contact  list,  the  new  entry  in  the 
contact  list  would  be  the  same  as  before,  except 

AA(K  +  2)  =  0.0 

AA(K  +  3)  =  0.0 

That  is,  the  normal  and  shear  contact  forces  would  be  set  to  zero. 

56.  Assuming  that  the  contact  between  elements  i  and  j  just 
discussed  was  the  first  contact  detected  for  element  i  ,  the  contact 
list  pointer  NAA(i)  is  set  equal  to  K  .  This  pointer  indicates  the 
position  in  memory  where  the  information  about  contacts  for  element  i 
begins.  Now  that  the  contact  data  for  this  contact  have  been  entered 
(and  the  entry  required  five  memory  locations),  the  next  available 
location  to  store  contact  list  data  is  determined  by  replacing  K  by 

K  +  5  (i.e.,  K  =  K  +  5).  Any  additional  detected  contacts  for  element  i 
are,  in  this  manner,  put  consecutively  into  the  contact  list;  each 
detected  contact  requiring  an  additional  five  memory  locations.  The 
index  K  is  incremented  by  five  following  each  entry.  After  all  con¬ 
tacts  for  element  i  have  been  detected  and  entered  into  the  contact 
list  (and  also  in  the  case  of  no  contacts  at  all  being  detected  for 
element  i) ,  the  end  of  this  contact  list  segment  (for  element  i)  is 
signified  by  embedding  a  zero  into  array  AA  at  location  K  +  1  ,  i.e., 

AA (K  +  1)  =  0.0 

The  embedded  zero  indicates  that  this  is  the  end  of  that  segment  of  the 
contact  list  pertaining  to  element  i  .  The  next  available  location  in 
array  AA  for  storing  contact  data  for  the  next  element  i  =  i  +  1  is 
now  located  at  K  =  K  +  2  .  Since  the  next  element  i  is  now  ready  to 
be  processed,  the  pointer  NA  ( i )  is  set  equal  to  K  .  The  just  described 
process  is  repeated  for  all  elements. 
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1 ni t ializat ion  of  contact  data  list 

57.  At  tho  outset  of  each  problem,  it  is  assumed  that  there  are 
no  contacts  between  the  various  disc  or  bar  elements.  As  mentioned  in 
the  preceding  paragraph,  zeros  that  are  embedded  in  the  contact  list  AA, 
at  the  location  of  the  value  of  the  element  pointer  NA(i)  ,  indicate  no 
contacts  for  element  i  .  Therefore,  at  the  beginning  of  a  problem 
containing  N  disc  and  bar  elements,  the  value  of  the  pointer  NA( i) 

is  set  to  i  yielding 

NA(1)  =  1 
NA  (2)  =  2 

NA(N)  =  N 

Since  a  zero  located  at  NA(i)  within  the  array  AA  indicates  no  con¬ 
tacts,  the  array  AA(K)  is  set  to  zero  (for  values  of  K  =  1,  N) .  In 
general,  the  first  time  cycle  in  MOTION  will  result  in  some  contacts. 

The  next  available  location  to  begin  storing  the  first  contact  within  AA 
is  at  N  +  1  .  The  new  pointer  value  for  the  first  contact  detected  will 
be  set  to  NA(i)  =  N  +  1  . 

58.  Therefore,  after  each  cycle  of  iteration  the  value  of  the 
pointer  NA(i)  is  the  location  within  the  array  AA  where  the  contact 
list  information  about  element  i  begins.  If  the  element  i  has  no 
contacts,  the  value  AA(NA(i))  will  be  zero,  signifying  the  end  of  the 
contact  list  segment  for  element  i  .  If  element  i  does  possess 
contacts,  the  value  AA(NA(i))  will  be  set  to  j  +  x^/10000  .  To  find 
the  end  of  the  contact  list  segment  for  element  i  ,  it  is  necessary  to 
skip  through  array  AA  by  increments  of  five  (starting  at  AA(NA(i))  until 
a  zero  is  found.  The  next  available  location  in  AA  will  possess  infor¬ 
mation  about  element  i  +  1  . 

Example  of 

c o ntac t  ii s t  inf ormat ion 

59.  Consider  the  simple  example  shown  in  Figure  7a.  At  the 
start  of  the  problem,  two  disc  elements  (2  and  1)  are  at  rest  on  bar 


element  1.  An  additional  disc  element  (5)  is  resting  upon  bar  element  4. 
Both  bar  elements  are  assumed  to  be  fixed,  i.e.,  they  are  prevented  from 
rotation  and  translation.  Assuming  that  the  only  forces  acting  on  this 
system  are  due  to  gravity,  the  subsequent  motions  of  the  disc  elements 
may  be  visualized  according  tu  the  sequence  indicated  by  Figure  7. 

Figure  8  shows  the  structure  of  the  contact  list  for  various  stages  of 
movement . 
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finis,  initially,  NA  (element  I)  =  1  ,  NA(2)  "  2  . . NA(5)  ~  1  ,  and 

the  contact  dal  a  list  (list  l,  figure  8)  contains  zeros  as  indicated  by 
the  symbol  K  (for  empty).  Sineo  there  are  live  elements  in  the  problem 
the  first  I  ive  storage  h  cal  ions  are  set  to  zero  (as  indicated  by  the 
t  ive  "F.'s"  in  list.  1).  The  next  available  location  for  rontart  storage 
(at  the  next,  or  in  this  ease  the  first  cycle)  is  location  6. 

62.  The  situation  for  the  1 irst  cycle  of  MOTION  is  shown  in 
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the  beginning  of  the  contact  information  for  contact  between  elements  1 
and  2.  Immediately  following  the  contact  information  for  contact  1-2  is 
the  information  for  contact  1-3  (starting  at  location  11  and  ending  at 
location  15).  The  zero  (symbol  E)  embedded  at  location  16  signifies 
that  there  are  no  more  contacts  for  element  1.  Contact  information  for 
element  2  begins  in  the  contact  list  at  location  17  (NA(2)  =  17). 
Location  17  contains  a  zero,  which  indicates  that  there  are  no  contacts 
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25  locations  have  been  used.  Tlius,  l  lie  next  available  locati  >n  foi¬ 
st  orage  is  26. 

64.  At  any  given  Lime  cycle,  it  is  necessary  to  retain  the 

contact  intormation  about  the  previous  time  cycle.  Only  by  retaining 
this  informal  ion  is  it  possible  to  check  whether  or  not  a  deLected 
contact  was  present  at  the  last  stage  or  il  the  contact  is  new.  (If  the 
contact  is  old,  the  normal  and  shear  forces  for  the  previous  cycle  are 
entered  into  the  new  contact  list;  otherwise,  zero  values  for  shear  and 
norma !  forces  are  entered.)  Figure  S  gives  the  new  contact  list  for 
cycle  2  (list  1).  List  1  is  exactly  the  same  as  1 ist  2  except  the  new 
contact  list  begins  at  location  26  and  ends  at  location  4r>.  1 1  is 

readily  apparent  that  all  of  these  "new"  contacts  entered  into  the  new 
contact  list  were  present  at  the  previous  cycle.  The  pointer  NA  for 

a  particular  element  is  not  updated  to  the  new  location  until  after  it 
is  used  to  check  whether  or  not  litis  same  contact  was  present  in  the 
previous  List.  For  example,  contact  1-2  has  been  identified  as  a  con¬ 
tact  for  cycle  2,  and  Lite  contact  information  is  to  be  entered  at 
locations  2b  through  30.  The  pointer  value  for  element  1.  N'A(l)  , 
still  has  a  value  of  6.  by  going  to  '.nation  6,  it  is  noted  that 
contact  1-2  was  present  In  the  previous  contact  list.  The  pointer  value 
N'A(l)  is  temporarily  updated  to  a  value  of  IL  so  as  to  be  directly 
pointing  to  the  next  location  in  the  previous  contact  list  where  contact 
information  about  element  1  was  stored.  Thus,  when  contact  1-3  is 
examined  (as  to  whether  or  not  it  was  in  the  previous  list),  the  pointer 
is  in  the  proper  position  for  quick  recognition  of  this  preexisting 
contact.  After  it  has  been  determined  that  all  contacts  for  element  1 
have  been  exhausted,  the  pointer  NA (11  is  set  to  the  new,  current 
value  of  26. 

65.  At  each  and  every  cycle  an  examination  is  completely  carried 
out  for  every  contact  given  in  the  poss ible  contact  list.  informal  ion 
stored  for  the  previous  cycle  is  not  used  at  all  for  determining  whether 
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or  not  a  contact  is  entered  into  the  new  contact  list'.*  Only  after  it 
has  been  decided  to  enter  a  contact  into  the  list  is  a  scanning  of  the 
previous  list  required  (and  then  only  to  decide  what  to  do  with  the 
shear  and  normal  forces).  The  method  for  updating  the  element  pointer 
NA(l)  simply  provides  an  efficient  method  for  checking  on  the  presence 
of  the  contact  in  question. 

66.  The  situation  for  cyclo  i  is  shown  in  list  4.  Again,  the 
organization  is  the  same  as  for  cycles  1  and  2  except  the  new  contact 
list  is  removed  hy  20  more  locations  (i.e.,  locations  46  through  65). 

The  situation  for  cycle  4  is  not  shown,  but  would  be  similar,  except  the 
new  contact  list  would  be  stored  in  locations  66  through  85.  This 
process  could  be  carried  out  indefinitely;  however,  at  some  point  memory 
would  become  exhausted.  Since  it  is  necessary  to  store  contact  infor¬ 
mation  only  for  the  previous  and  present  cycles,  a  scheme  was  developed 
for  recovering  all  available  memory. 


There  is  an  option  available  that  precludes  the  searching  for  con¬ 
tacts  at  each  and  every  cycle.  By  entering  a  suitable  variable 
(Command  I  in  the  input  Phase  (Figure  9)),  it  is  possible  to  stipulate 
that  complete  searches  for  contacts  will  be  mad e  only  after  every  so 
many  eycLes.  Unless  otherwise  stipulated  through  these  options,  a 
complete  search  is  made  for  those  element  contacts  in  the  possib le 
contact  list  at  each  cycle.  If  the  option  is  used,  the  program  as¬ 
sumes  that  the  current  real  contact  list  is  to  be  used  over  and  over; 
the  shear  and  normal  forces  are  being  constantly  updated.  During  this 
time,  the  contact  coordinates  and  contact  plane  angles  are  properly 
computed.  The  use  of  this  option  can  be  dangerous  since  new  contacts 
wilL  not  bo  detected  until  a  complete  search  is  called  for.  A  com¬ 
plete  search  is  always  made  upon  each  entry  into  MOTION.  The  use  of 
this  option  can  significantly  reduce  computer  processing  t  ime  since 
it  negates  the  need  for  contact  searches  and  for  creating  a  new  con¬ 
tact  list. 
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t)  7 .  Suppose  that  12U  storage  locations*  have  been  reserved  for 


contact  information  as  indicated  in  Figure  8.  While  processing  in 
cycle  1,  iL  is  apparent  that  locations  1  Lhrough  25  are  needed  for  the. 
present  and  the  previous  contact  lists.  Thus,  locations  26  Lhrough  120 
are  available  as  "free"  memory  for  future  cycles.  While  in  cycle  2,  the 
required  memory  locations  are  6  through  45  (old  plus  new  contact  list). 
Therefore,  the  free  memory  consists  of  locations  46  through  120  plus  1 
through  5.  At  cycle  3,  the  free  memory  consist's  ol  locations  66  through 
120  plus  1  through  25.  During  cycle  4,  the  free  memory  extends  from 
locations  86  through  120  plus  1  throng))  45.  Al  cycle  5,  the  free 
memory  consists  of  locations  106  through  120  plus  1  through  65,  Thus, 
as  cycle  6  is  entered  (and  will  require  20  locations  to  store  the  new 
contact  list),  iL  is  apparent  that  the  reserved  memory  space  of  1.20  will 


The  size  of  common  memory  for  the  variables  used  in  program  DISC  is 
set  through  the  use  of  a  Fortran  PARAMETER  statement.  As  seen  from  an 
examination  of  the  program  listing,  the  space  allocated  for  the  contact 
list  AA  is  set  to  M  x  24  ,  where  M  is  the  maximum  number  of  elements 
permitted  for  a  problem.  In  practice,  M  will  often  be  greater  than 
the  actual  number  of  elements  being  considered,  thus  extra  memory  space 
is  generally  allocated.  The  common  dimensions  of  M  x  24  should  be 
adequate  for  most  problems.  For  example,  consider  a  collection  of 
closely  packed  uniform  disc  elements,  x  discs  wide  by  y  discs  high. 
The  total  number  of  elements  is  xy  .  The  maximum  total  number  of  con¬ 
tacts  for  this  collection  may  be  shown  to  be  )xy  -  2(x  +  v)  +  1  • 

Since  five  storage  spaces  are  required  for  each  current  contact  (plus 
one  additional  space  for  each  contact  to  store  the  zero  or  "F."  marks) 
the  storage  S  required  to  contain  the  current  contac*  list  is 

S  =  5[ (3xy  -  2(x  +  v)  +  1 ) |  +  xv 
=  16xv  -  1.0 (x  +  y)  +5 

Neglecting  the  last  two  terms,  S  equals  .16  times  the  total  number  of 
elements.  Technically,  twice  this  amount  of  storage'  is  required  at  any 
given  cycle  since  the  last  contact  list  must  be  scanned  to  prepare  the 
current:  contact  list.  In  actuality,  twice  the  amount  is  not  needed 
because  as  soon  as  the  new  contacts  are  entered  into  the  » urrent  con¬ 
tact  list,  there  is  no  longer  any  need  for  the  corresponding,  previous 
contact  (as  long  as  the  sequence  of  insertion  remains  the  satin  -i  •  At 
present,  it  is  felt  that  a  reserved  space  of  M  24  should  he  ode- 
quale  for  most  problems. 


I'll  is  condition 


be  exceeded  unless  the  I  rent  part  et  memory  is  reused, 
can  be  avoided  by  romput  in;',  the  amount  of  tree  memory  still  available  up 
to  location  120  niter  each  entry  into  the  contact  list.  Since  the  next 
entry  could  require  five  additional  storage  locations,  the  beginning  of 
the  next  entry  is  directed  to  location  1  if  sufficient  memory  does  not 
exist.  Therefore,  as  shown  tor  cycle  b,  the  first  detected  contact  1-2 
is  placed  in  locations  106  through  110,  and  the  next  contact  placed  in 
locations  111  through  116.  The  five  storage  locations  available  before 
location  120  are  sufficient  for  storage  of  one  more  contact.  However, 
contact  1-3  was  the  last  contact  for  element  1.  Therefore,  a  zero  ("E") 
is  placed  at  Location  1L6,  signifying  the  end  of  contacts  for  element  1. 
There  are  now  only  four  available  locations  before  the  end  of  the  re¬ 
served  memory  and,  therefore,  not  enough  locations  to  store  information 
for  a  detected  contact.  The  program  now  directs  that  storage  begin  at 
location  1.  As  shown  for  cycle  6,  a  zero  ("E")  is  placed  in  location  1 
(signifying  that  element  2  has  no,  as  yet,  undetected  contacts).  The 
pointer  NA(2)  is  set  to  the  value  of  1.  Locations  117  through  120  are 
not  used. 

68.  As  long  as  the  contacts  do  not  change,  the  calculations  will 
proceed  as  described  In  cycles  1  through  6,  each  cycle  requiring  an 
additional  20  new  storage  locations  of  free  memory.  Cycle  457  is  the 
situation  just  before  element  2  rolls  off  bar  element  1.  Figure  8 
presents  the  possible  contact  list  existing  for  the  conditions  shown  in 
Figure  7b.  The  interpretation  of  this  list  is  as  follows;  elements  2 
and  3  tmay  contact  elem<  nt  1,  element  4  may  contact  element  2,  element  3 
has  no  undetected  contacts,  element  5  may  contact  element  4,  and  element 
has  no  new  contacts.  The  possible  contacts  1-2  and  2-4  are  included  in 
the  possible  contact  list;  however,  they  will  be  excluded  as  real  con¬ 
tacts  in  MOTION  since  the  situation  shown  in  Figure  7b  is  for  the  condit 
where  element  2  is  just  breaking  contact  with  element  1  and  bet  ore 
contact  with  element  4.  Figure  8  also  shows  the  new  contact  list  for 
cycle  458.  Note  that  there  are  now  only  two  contacts  (contacts  1-3  ■  d 
4-5).  The  contact  list  for  cycle  459  is  just  as  it  was  for  cvcle  43, s 
except  the  new  contact  is  stored  an  additional  15  locations  turtlur 
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along  Ln  memory.  Cycle  679  is  the  situation  just  before  element  2  hits 
element  4.  ln  addition.  Figure  8  presents  the  possible  contact  list  and 
the  resulting  new  contact  lists  (cycles  680,  681,  and  1085).  Notice 
that  again  there  are  three  contacts  (contacts  1-3,  2-4,  and  4-5). 

69.  In  Figure  8,  the  possible  contact  list  and  the  contact  data 

list  are  given  for  some  much  later  time  (say  cycle  2002).  These  lists 

are  for  the  conditions  depicted  in  Figure  7d.  By  this  time,  all  of  the 

disc  elements  are  in  contact  and  proceeding  to  roll  off  the  end  of 

element  4.  At  this  stage,  five  contacts  are  detected  (contacts  2-3, 

2-4,  2-5,  3-4,  and  4-5).  Thus,  30  additional  storage  locations  are 

required  to  store  the  contacts  at  each  cycle. 

Calculations  of 
forces  on  t he  elements 

70.  As  each  element  is,  in  turn,  examined  for  contacts,  the  sum 
of  all  forces  acting  on  the  element  is  computed.  Preceding  each  time 
cycle  of  iteration  the  force  sum  is  set  equal  to  the  applied  plus  gravity 
loads,  i . e .  , 

F*  =  vertical  applied  load  minus  the  gravity  load 
ysum  n  j 

F1  =  horizontal  applied  load 
xsum 

Ml  =  applied  moment 
sum 

In  the  event  that  the  element  being  examined  possesses  no  contacts,  the 
computation  proceeds  directly  to  the  calculation  of  accelerations  and 
displacements  (see  Equations  6  through  9,  Part  1 L) .  if  contacts  are 
found  for  the  element  being  considered,  then  the  force  sums  are  modified 
in  accordance  with  Equations  1  through  5.  The  force  sums  are  updated 
for  both  elements  in  contact.  For  example,  suppose  a  situation  exists 
in  which  element  3  makes  contact  with  elements  5  and  7.  When  element  3 
is  examined  for  its  contacts,  the  force  sums  will  be  updated  twice  (once 
for  element  5  and  again  for  element  7).  During  this  process,  the  force 
sums  for  elements  5  and  7  are  also  updated.  Then  when  element  5  is 
examined  and  it  Is  determined  that  it  has  no,  as  yet,  undetected  contacts, 


the  computation  will  proceed  directly  to  Equation  6  since  the  force  sum 
for  element  5  was  updated  at  the  time  when  element  3  was  examined.  The 
same  holds  for  element  7. 


71.  The  computations  of  liquations  8  and  9  are  particularly 
important.  Equation  8  yields  the  new  incremental  element  centroid 
displacements,  and  Equation  9  yields  the  new  coordinates  of  the  element 
centroid.  It  is  these  quantities  that  are  required  by  Equation  1  for 
the  next  iteration  cycle. 

Other  Subroutines 

Su  b  routine  PLTFOR 

72.  Subroutine  PLTFOR  is  used  to  output  the  forces  acting  at  the 
contacts  between  elements.  Recalling  that  the  contact  list  is  struc¬ 
tured  such  that  each  element  has  at  least  one  entry  (if  no  contacts  were 
detected)  and  further  that  the  storage  is  in  sequential  order,  it  is 
possible  to,  at  any  stage,  thread  through  the  contact  list  to  obtain  the 
contact  Lorces.  For  each  element  i  that  possesses  contacts,  the 
contact  list  contains: 

a.  The  number  j  of  the  element  being  contacted  by  i  . 

b.  The  ID  number  for  the  contact  (two  bar  elements  in 
contact  may  have  up  to  four  contact  points). 

c.  The  contact  point  coordinates,  x  and  y 

c  c 

d.  The  normal  and  shear  forces  acting  at  the  contact,  F 

and  F^  (stored  as  AA(K  +  2)  and  AA(K  +  3)  in  memory). 

e.  The  angle  to  the  contact  plane,  B  . 

73.  These  forces  are  output  graphically.  A  pass  is  first  made 
through  the  contact  list  to  determine  the  largest  absolute  value  of  the 
shear  or  normal  force  in  the  list.  This  quantity  is  used  to  calculate  a 
scale  factor,  and  the  vector  representing  this  force  will  be  20  screen 
units  long.  The  other  vectors  drawn  to  represent  the  shear  and  normal 
forces  will  be  proportional.  Thus,  after  each  call  to  PLTFOR,  a  graphical 
representation  of  the  contact  forces  are  drawn  on  the  screen.  The 
contact  coordinates  X(  and  v  are  used  to  locate  the  position  of  the 
vectors  to  be  drawn,  the  shear  force  is  drawn  at  the  angle  B  ,  and  the 
normal  force  perpendicular  to  the  shear  force.  The  user  also  has  the 
option  of  obtain  in;;  a  printout  of  the  contact  forces  (Figure  9  or  10). 


Sub rou Line  LOOP 

74.  Subroutine  LOOP  is  used  to  check  whether  or  not  the  cross¬ 
hair  cursor  is  located  in  the  vicinity  of  an  element  centroid.  In  some 
cases,  this  test  is  used  to  modify  certain  commands. 

Su broutine  INTERVAL 

75.  Subroutine  INTERVAL  is  used  to  change  Lhe  frequency  of 
returning  control  to  the  terminal.  The  normal  method  is  to  perform 
250  time  steps  of  calculations  before  halting  calculations.  At  each 
halt  the  user  directs  what  action  is  to  be  taken  next. 

Subjrou_t_i_ne  GRIP 

76.  Subroutine  GRID  causes  a  grid  to  be  drawn  on  the  screen. 
Subroutine  SAVE 

77.  Subroutine.  SAVE  may  be  called  at  any  time  by  OUTPUT  to  store 
the  current  element  geometry  in  memory.  Subroutine  SAVE  is  automatically 
called  on  each  exit  from  INPUT.  This  feature  provides  the  user  a  means 

of  reworking  the  same  geometry  with,  perhaps,  a  different  set  of  parameters. 
S ubrou tine  CIRCLE 

78.  Subroutine  CIRCLE  aids  in  the  graphical  output  of  d  isc  and 
bar  elements. 

Sub rot i t in e  SCREEN 

79.  Subroutine  SCREEN  is  used  to  position  the  alpha  cursor  on 

the  screen  following  input  queries  and  some  output  functions.  It  prevents 
overwriting  certain  areas  of  the  screen. 

Su  broutine  GEN 

80.  This  subroutine  may  be  called  during  the  INPUT  phase.  Its 
function  is  to  create  a  collection  of  disc  elements  touching  one  another 
in  a  dense  packing. 

Subroutine  VECTOR 

81.  This  subroutine  is  used  to  vary  the  form  of  Lhe  graphical 
output.  The  norma  1  mode  of  display  is  to  draw  the  elements  on  tin.1 
screen  at  their  current  positions.  Subroutine  VECTOR  causes  the  display 
to  track  Lhe  element  centroids. 


Sub  rout  ino  WF.  I  pHT 

82.  Subroutine  WEICHT  may  be  used  to  modify  the  weight  of  any  or 
all  elements. 

Subrou t  iaie_  DAMP 

83.  Subroutine  DAMP  is  used  to  change  the  damping  constant  K  . 


I’ APT  IV: 


I’|««;i<a:i  ni'i.it.Yi  m:. 


■84  .  Ihe  eomput  or  program  DISC  w ,  i wi  it  ten  in  -.u.  h  .1  1  .isli  i  mi  .1;. 
to  keep  the  data  prep, amt  ion  simple  .uni  minimal.  A:.  1  result  ,  n:u  >t 
computer  graphics  was  emp  1  oved  to  a  large  extent.  Mi > • . t  ot  tin.-  cnmi'uml: 
lU'i-i'ssarv  to  opi-ialt-  DISC  consist  ol  single  ktV-.lioi-.es  .1  tin.  Tek  t  t  011  i 
terminal.  Figures  9  ami  ID  arc  the  computet  iim-i  gui, tot  tie 
operation  ol  DISC.  Figure  9  gives  a  more  Detailed  d  i  •>  us-  : on  -1  1  lie 

user  commands;  Figure  !()  is  an  abbreviated  command  I  ist  ini'.  '.-.'lien  one 
I  irst  uses  the  program  DISC,  the  inlormat  ion  given  in  Figure  '*  should 
guide  the  user  in  tin-  operation  ol  the  program.  Alter  the  user  become 
familiar  with  the  program,  the  abbreviated  command  li  t  (Figure  1  <  1  > 
sliiui  Id  su I  f  ice. 

Sr).  According  to  the  user's  guides,  all  input  is  initiated  by 
entering  a  single  character,  following  the  appearance  ol  the  crosshair 
cursor.  The  crosshair  cursor  consists  of  two  perpend  icul ar  lines  cros 
the  screen  ol  tile  terminal.  The  position  ol  the  inlei  sect.  ions  ot  the 
two  Lilli’s  is  controlled  by  a  set  ol  thumbwheel.-,  on  tin  terminal.  In 
many  instances,  the  entry  ol  the  single  character  is  sufficient  to 
accomplish  whatever  is  desired.  In  other  instances,  certain  parameter 
(or  numbers)  associated  with  the  single  character  command  are  required 
if  tile  command  requires  one  or  pore  parameters,  the  alpha  cursor  (a 
small  matrix  of  dots)  will  appear  on  the  screen.  As  soon  as  tile  alpha 
cursor  appears,  the  input  parameters  may  he  typed. 

8h .  Tin’  detailed  user's  guide  (Figure  9)  is  complete  regarding 
the  meaning  of  the  various  commands  and,  as  such,  will  not  he  discusse 
further.  Three  examples  will,  instead,  lie  described  in  order  to  illus 
trate  the  program. 

87 .  Figure  II  is  a  series  ol  computer  drawings  that  describe  a 
very  simple  problem.  As  soon  as  program  DISC,  is  called  into  exeeul  ion 
the  message  "INPUT  P1IASK"  is  written  on  the  screen  <  Figure  11a).  The 
col leet ion  of  throe  elements  was  created  in  t he  following  order: 
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174  t  £  Ends  compute  phase.  Calls  input  phase 
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ImPUT  phase 


COMPUTE  PHASE 


(Shfi't  4  of  10) 
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this  time  mi,  it  is  not  possible  to  cause  the  small  disc 
separate-  from  the  bar  element.  file  lorees  (in  both  normal  ami  shear 
ilirei  t  ions)  are  ontputteil  to  the  upper  r  ight  corner  ol  the  screen  b\  tin 
comma  nil  "0"  I  ol  lowed  by  "II."  Thus,  lor  the  contact  between  elements  1 
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hi.  At  this  point,  the  bar  element  is  released  from  its  i  ixitv 
in  all  directions  through  the  command  "U"  (with  the  crosshairs  centered 
on  the  bar  element).  This  command  is  followed  hi’  restraining  the  bar 
element  from  movement  in  both  the  horizontal  and  vertical  directions  be 
striking  "X"  I  ol  lowed  by  "V"  (again  with  Lite  erosshair  cursor  on  c  I  emenl 
The  bar  element  is  now  prevented  from  translating  but  is  free  to  rotate 
about  its  center  of  giavitv.  fo  attain  begin  culcnlal  ions  ot  motion,  tin 
key  is  struck.  him  e  tin-  large  disc  exert  s  a  counterc  lockwise 
moment  on  the  bar,  the  sv  .ran  begins  to  rotate  as  shown  in  Figure  lie. 

hi’.  Repealed  "<!"  oBiiimd;;  lead  to  the  situation  shown  in  Figure 
and  Figure  lig.  The  bar  e  on  t  i  nm  s  to  t  ot  ate,  am!  the  d  isc  elements  rol 
down  the  p'nne.  Ihe  lines  joininr,  the  disc  eLemcnt  centers  to  Lite 
periphery  ol  the  discs  indicate  the  amount  of  rotation  experienced  be 
the  discs  (Figure  Mg,).  At  a  later  time  (Figure  I  ill),  the  large!  disc  . 
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shows  the  resulting  mol  ions  alter  striking  key  "(1 . "  Notice  that  disc  h 
proceeds  downward  to  the  left.  Disc  7  starts  upward  to  the  right, 
strikes  Lite  upper  boundary,  and  turns  downward  to  the  right.  Following 
the  first  calculation  cycle  of  250  time  steps,  output  is  converted  to 
the  "vector"  mode  by  hitting  key  When  in  the  vector  mode,  the 

elements  are  no  longer  drawn.  1'he  output  now  consists  only  of  the 
movement  of  the  element  centers.  Klement  7  proceeds  to  rebound  off  the 
lower  boundary,  then  strikes  the  lower  end  of  element  2,  and  exits 
through  the  opening  (Figure  1 J b ) .  Element  6,  however,  does  not  quickly 
find  an  exit  through  the  opening.  In  Figure  13c,  element  6  suffers  a 

total  of  21  rebounds  before  hitting  the  rounded  portion  of  element  2. 
Element  6  then  is  deflected  in  an  almost  horizontal  direction  to  the 
left  wall  and  shoots  across  and  out  through  the  opening  (almost  dead 
center) . 

96.  The  purpose  of  the  three  example  problems  is  to  demonstrate 
the  flexibility  and  range  of  this  program.  Problems  involving  quasi¬ 
static  behavior  and  large  rapid  motions  are  included  within  the  same 
framework.  Indeed,  at  one  time,  this  program  contained  options  for 
including  "inverse  square  of  the  distance"  forces  so  that  the  user  could 
compute  orbits  of  gravitational  bodies.  These  options  are  not  germane 
to  the  general  fields  of  geotechnical  engineering  and,  as  such,  are  not 
present  in  this  version.  However,  the  program  DISC  is  basically  a  tool 
that  will  solve  many  fundamental  problems  in  rigid  body  mechanics. 
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II  sli.il  1  remain  tor  others  to  extend  the  dist  ini  t  <  lenient 


method  to  ;::ut  e  meaningful  situations.  It  is  the  writer's  opinion  that 
this  t  vpe  of  analvsis  will  eventual  lv  lie  adopted  in  many  areas  ol  1 1  i  - 
luerini;.  i  In  ooneept  of  analyzing  discrete  partieles  in  an  explicit 
scheme  is  attr.u  t  ive  to  the  formulations  of  goot  ee  hn  i  ea  I  prolilems. 

Reeoiumendat  ions 

lt)t).  It  must  hi'  emphasized  that  the  computer  program  DISC  was 
writ  tint  to  be  quite  general.  it  a  subsequent  user  of  this  program  has 
to  analyze  a  particular  problem,  lie  should  make  an  attempt  to  modify  the 
basic,  general  program  to  Itis  needs.  For  example,  if  the  problem  re¬ 
quires  only  a  quasi-static  analysis,  many  savings  could  be  accomplished 
hv  fewer  searches  for  near  contacts.  Researchers  generally  tend  to 
direct  new  concepts  to  a  large  audience  and  always  try  to  account  I  or 
anv  possible  upplicat  inn  of  their  developments,  A  user  of  this  research 
would  he  well  advised  to  simply  store  a  copy  of  the  original  program 
DISC  within  his  files  and  then  mod i 1 y  and  enhance  whichever  portions  ol 
DISC  that  fit  his  app 1 icat ions .  That  is,  adapt  the  program  to  the 
pa  r  t  i  i  n  1  ar  need  . 
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